Molecular electronics and first-principles methods 
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We discuss the key steps that have to be followed to calculate coherent quantum transport in 
molecular and atomic-scale systems, making emphasis on the ab-initio Gaussian Embedded Cluster 
Method recently developed by the authors. We present various results on a simple system such as 
a clean Au nanocontact and the same nanocontact in the presence of hydrogen that illustrate the 
applicability of this method in the study and interpretation of a large range of experiments in the 
field of molecular electronics. 



I. INTRODUCTION 

Nanoscale electronics constitutes the backbone of 
present and future technological advances or, in other 
words, of Nanotechnology. Ultimately, the functionality 
of electronic devices will rely on the conduction proper- 
ties of nanoscopic regions composed of a number of atoms 
that can range typically from several thousands down to 
a single one. Probably, the most promising research field 
within nanoscale electronics is what is known as molecu- 
lar electronics. The main idea behind molecular electron- 
ics is the possibility that functional units can be built out 
of very stable and well-characterized molecules 1 - such as 
fullerenes 2 or carbon nanotubes 3 ^. The simplest molecu- 
lar device consists of two large metallic electrodes, several 
nanometers apart, joined by a molecule or molecules an- 
chored to them-. We will name this system molecular 
bridge. Alternatively, the electrodes can be simply con- 
nected by a chain of atoms of the same element as the 
electrodes&t These bridges, generically known as atomic 
contacts or metallic nanocontacts, are not expected to be 
of any practical application, but constitute an excellent 
benchmark for us to learn about the world of electrical 
transport at the atomic scale. 

The design and fabrication of electronic devices at the 
molecular and atomic scale poses new challenges for the- 
orists who must develop and implement new techniques 
to address the upcoming problems. The basics to cal- 
culate the zero-bias, zero-temperature conductance G of 
a molecular bridge or metallic nanocontact were estab- 
lished by Landauer long before the concept of molecular 
electronics was commonplace. In Landauer 's formalism 
G is simply related to the quantum mechanical transmis- 
sion probability T of electrons at the Fermi level to go 
from one electrode to the other— (we assume spin degen- 
eracy): 

G = G Q T(E F ) = \t{E ¥ ). (1) 
The transmission can be easily estimated on generic con- 



siderations for metallic nanocontact o^i 1 ? , but it is much 
more difficult to do so for molecular bridges. The rea- 
son is simple: We do not know a priori where the elec- 
trode Fermi energy E-p lies with respect to the molecular 
levels. The positioning of Ep depends on two closely- 
related factors: (i) the coupling or hybridization of the 
molecular levels with the free-electron levels in the elec- 
trodes and (ii) the charge transfer between molecule and 
electrodes. The conductance is thus strongly dependent 
on the particular molecule, the detailed atomic arrange- 
ment of the electrodes where the molecule binds, and the 
chemical nature of the various elements at play. In or- 
der to give an answer to this problem we have to rely 
on first-principles or ab-initio calculations. To calculate 
the atomic arrangement of the electrodes or the way the 
molecule binds to the electrodes is a major problem in 
itself that deserves a whole study on its own and will not 
be given further consideration in these notes. Still, even 
if we ignore this important issue, to implement Landauer 
formalism requires to know the electronic structure of a 
finite region embedded in an infinite system with no par- 
ticular symmetry or periodicity. This highly non-trivial 
problem is what we discuss in what follows. 



II. THE BASICS 

The main advantage of the numerical implementa- 
tion that we discuss here with respect to similar pro- 
posals that have appeared in the literatureiiiiSii 3 * 1 ^ is 
the use of a standard quan tum chem istry code such as 
GAUSSIAN98^ (see Refs. 1 1 6l 1 7l 1 8lfT^I . The GAUS- 
SIAN98 code provides a versatile method to perfom first- 
principles or ab-initio calculations of clusters, incorporat- 
ing the major advancements in the field in terms of func- 
tionals, basis sets, pseudopotentials, etc.. The procedure 
goes as follows. A standard electronic structure calcula- 
tion of the region that includes the molecule and a signif- 
icant part of the electrodes is performed [see Fig. da)]. 
This calculation can be performed within any mean-field- 
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like approach: Hartree-Fock or density functional (DF) 
theory in any of its multiple approximations. As far 
as transport is concerned, the self-consistent hamilto- 
nian H (or Fock matrix F) of this central cluster or 
supermolecule contains the relevant information. The 
retarded(advanced) Green's functions associated with F 
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(2) 



needs to incorporate the rest of the infinite electrodes: 
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In this expression 
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where Er(E l ) denotes a self-energy operator that ac- 
counts for the part of the right (left) semi-infinite elec- 
trode which has not been included in the calculation, 
and J is the unity matrix. In a non-orthogonal basis the 
Green function takes the form 
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where S is the overlap matrix. To simplify the notation 
we have dropped the energy argument in the Green func- 
tion and the selfenergy. The added self-energy matrices 
can only be explicitly calculated in ideal situations. We 
choose to describe the bulk electrode with a Bethe lat- 
tice tight-binding model with the coordination and pa- 
rameters appropriate for the electrodes included in the 
cluster. Figure^c) depicts schematically a Bethe lattice 
of coordination three in two dimensions. The advantage 
of choosing a Bethe lattice resides in that it reproduces 
fairly well the bulk density of states of most commonly 
used metallic electrodes. For each atom in the outer 
planes of the cluster, we choose to add a branch of the 
Bethe lattice in the direction of any missing bulk atom 
(including those missing in the same plane). In Fig.^b) 
the directions in which Caylee tree branches are added 
are indicated by bigger atoms which represent the first 
atom of the branch in that direction. Assuming that the 
most important structural details of the electrode are in- 
cluded in the central cluster, the Bethe lattices should 
have no other relevance than that of introducing a generic 
reservoir. 

The above definition of the Green function is not a 
standard one. As many other ab-initio implementations, 
GAUSSIAN98 makes use of non-orthogonal localized or- 
bital basis sets. The appearance of S in the definition 
of the Green function guarantees that the total charge 
of the system N can be obtained through the standard 
expression used in quantum chemistry: 



N = Tt[P-S], 



(6) 



where the trace runs over all the atomic orbitals of the 
cluster and where P is the density matrix defined in a 




FIG. 1: (a) A molecular bridge model where a short capped 
nanotube is contacted by two (001) parallel surfaces. The 
standard self-consistent procedure performed initially by 
GAUSSIAN98 is also shown schematically on the right, (b) 
The same system where the first phantom atoms of the Bethe 
lattices are also depicted. The self-consistent procedure when 
the Bethe lattices are included is also shown schematically 
on top. (c) A finite section of a Bethe lattice model in two 
dimensions with coordination three. 



standard way in terms of the retarded Green function as 



P = 

7T 



Im 
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E-p is the Fermi energy which can be set by imposing 
overall charge neutrality in the cluster. The integral in 
Eq. J7J can be efficiently calculated along a contour in 
the complex planeiii^S. As explained in the next sec- 
tion, we force GAUSSIAN98 to evaluate F using the den- 
sity matrix defined in Eq. instead of the standard one 
obtained from filling up the molecular orbitals obtained 
from diagonalizing the Fock matrbiSi. It is also inter- 
esting to note that the applicability of this approach, 
which we name the Gaussian Embedded Cluster Method 
(GECM) can be easily implemented in any quantum- 
chemistry or ab-initio code as long as it is based on a 
localized orbital expansion of the self-consistent, single- 
electron wavef unctions. 

The transmission probability that appears in Eq. ^ 
can be calculated through the general expression^ 



T{E) =Tr[f L G (+) f R G ( - ) ], 



(8) 



where Tr denotes the trace over all the orbitals of the 
cluster and where the operators Ir(l) are given by 
i(£r(l) — Sj^L-j). Finally, in order to single out the con- 
tribution of individual channels to the current, one can 
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FIG. 2: Schematic picture of the development of a potential 
drop between electrodes separated by a barrier on imposing an 
electrochemical potential difference between them, (a) Before 
self-consistency begins the chemical potential in electrodes is 
different by an amount eV and the bottom of the conduc- 
tion bands is the same for both electrodes. This situation 
corresponds to a diffusion problem, (b) After self-consistency 
has been achieved the chemical potential difference practi- 
cally vanishes in favor of an electrostatic potential difference 
between eletrodes. 



diagonalize the transmission matrix. While the size of 
the matrix [f L G^t R G^\ can be as large as desired, 
the number of eigenvalues with a significant contribution 
will be typically much smaller, being determined by the 
narrowest part of the nanocontact or the molecule. The 
symmetry of each channel can be identified by looking at 
its weight on the atomic orbitals of the cluster. 

Notice that, although the above expression provide us 
with the transmission as a function of energy, only the 
value of T at the Fermi energy, T(Ep), is strictly cor- 
rect (within the approximations inherent to DF theory). 
Away from Ep, T(E) is only indicative. For instance, if 
a gate voltage is applied and the system gets charged or 
discharged, T needs to be recalculated for the new Fermi 
energy. To obtain the current I at non-zero temperature 
and finite bias voltage V one integrates the transmission 
probability: 

2e f°° 

I = J T(E, V)[f L (E - eV/2) - f R (E + eV/2)]dE, 

where /l and /r are the Fermi distribution functions of 
the left and right electrodes, respectively. Notice, once 
again, that this innocent-looking expression is not a triv- 
ial generalization of Eq. ^extended to finite bias voltages 
and temperature. The transmission coefficient that ap- 
pears in the integrand T{E, V) depends on temperature, 
but, most importantly, depends on V . To calculate this 
transmission coefficient, knowledge of the Fock matrix in 



the presence of a bias voltage is required which, in turn, 
depends on the density matrix out of equilibrium. Based 
on standard non-equilibrium Green function techniques, 
it has been repeatedly shown in the literature 8 that 



r = -— i dE [s- 1 g < s- 1 ] 



(10) 



where 



G<(E) = iG^-\E) \f h (E)t h (E) + f R (E)T R (E) 



G (+) (£). 
(11) 

There are several technical and conceptual issues regard- 
ing the computation of the above equations that need 
some discussion. We refer the reader to the literature for 
a detailed account of most of themi^iiiii. Here we will 
simply mention two of them. Maybe the most important 
conceptual issue is that related to the fact that using DF 
theory for non-equilibrium situations is not justified: DF 
theory is a ground-state theory. We will make use of it 
simply because we do not know a better way to deal with 
out-of-equilibrium problems in an operative way. Second, 
there is technical issue that, in our view, has not received 
due attention in the literature and that we find it worth 
a detailed discussion^: The out-of-equilibrium electro- 
statics. In most self-consistent approachesii*i^*i^iiS one 
solves the Poisson equation with boundary conditions ap- 
propriate for the electrode geometries. This is equivalent 
to imposing an external electric field and it can only 
be done for simple geometries such as infinite planes. 
The Keldysh formalism, however, does not necessarily 
require to deal with the Poisson equation if a significant 
part of the metallic electrodes has already been included 
in the central cluster. Thus, one can consider realistic 
electrode geometries if necessary. One simply imposes 
an electrochemical potential difference between electrodes 
Mec — Mec = e ^ which is actually what a battery does. 
There are two contributions to fj, ec : 

Mec — Mchcmical T" /^electrostatic i 

where we define 



(12) 



Mchcmical — 

lim SE COTe /5N 



/^electrostatic 



SN^O 

lim 5E ee /SN. 

SN^O 



E COTe is the standard energy contribution to the total en- 
ergy from the cores of the atoms and E ee is the contribu- 
tion from the electron-electron interaction^. Specifically 

SE COTC = tr[5P ■ h core ] 
6E CC = tr[5P-h ee ], 



where h c 



and h ee are the core and interaction ma- 



trices composing the Fock matrix. In the selfconsistent 
procedure one electrode gets charged and the other dis- 
charged (by the same amount in case of mirror symme- 
try). Once self-consistency has been achieved, an elec- 
trostatic potential difference V between atoms one or 
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FIG. 3: Average on-site energies for all the atoms in the Au 
nanocontact shown in the inset. White dots have been ob- 
tained at zero bias and black dots at 5 V. 



FIG. 4: Transmission versus energy for 1^=0 (continuous line) 
and 5 V (dashed line) for the Au nanocontact shown in Figure 
3. 



two layers inside opposite electrodes develops while their 
chemical potentials are similar, i.e., they remain neutral. 
On the other hand, the electrostatic potential drop will 
be smaller than V between atoms at the surface of op- 
posite electrodes since they carry the charges. We show 
this process schematically in Fig. [5] 

III. EXAMPLE: AU NANOCONTACTS 

A. Clean nanocontacts 

Calculations that illustrate the methodology described 
above have been carried out for the Au nanocontact 
shown in the inset of Fig. Two (111) planes containing 
19 atoms each plus a four-atom chain. Interatomic bulk 
distances have been taken for the whole cluster (2.88A), 
although one should keep in mind that an analysis of the 
structural stability of the nanocontact model is an im- 
portant issue when one is interested in the interpretation 
of experimental data. The results shown next correspond 
to zero temperature and for the DF calculations we have 
used the Becke's three-parameter hybrid functional using 
the Lee, Yang and Parr correlation functional (B3LYP) 23 
together with the semilocal shape consistent pseudopo- 
tential (SCPP) and minimal basis sets of Christiansen 
et oZi 24 i 25 . Fig. |31 shows the average on-site energies of 
the 5d6s6p orbitals for all the atoms of the nanocontact 
at zero and 5 V bias. This magnitude reflects the local 
self-consistent electrostatic potential on each atom. The 
results represented by white dots correspond to zero bias 
where one can see the characteristic potential profile due 
to the parallel surfaces reflected on the chain atoms ly- 
ing across the vacuum between surfaces. The black dots 
correspond to a 5 V bias. One can see how an electro- 
static potential drop of almost 5 V has developed between 
plates without having imposed an external field, just an 



electrochemical potential difference. The major drop in 
the electrostatic potential occurs at the chain-plane con- 
tacts. In the case of zero bias the results are symmetric 
with respect to the geometric center of symmetry, as ex- 
pected. Instead, a similar symmetry is absent for 5 V. 
Namely, whereas the potential drop between the left elec- 
trode and the first atom in the chain is 1.92 eV, it is only 
of 0.85 eV between the chain end and the right electrode. 
This feature seems common to all previously reported 
results 13 '—- and reflects bulk band structure features. The 
transmission T for zero and 5 V is shown in Fig^ T oscil- 
lates around the Fermi energy with an upper limit value 
of one for both cases. This limit is well documented the- 
oretical and experimentally for Au nanocontacts^ since 
there is only one s-like channel around the Fermi energy 
for Au chains. The oscillations are due to scattering at 
the electrode-chain contact. This explains why conduc- 
tance histograms exhibit the lowest peak slightly below 
G = 2e 2 /h 2S . It is worth noting that although the total 
current calculated by integrating T(E, V = 0) in a win- 
dow [-V/2, V/2] does not differ much from that obtained 
with the full non-equilibrium approach, the differential 
conductance is significantly different. For instance, the 
gap below —0.5 eV is partially filled at finite bias. 



B. H2 molecular bridges in Au nanocontacts 

We finally examine the possible consequences of adding 
hydrogen (H2) to a Au nanocontact. A similar experi- 
mental analysis has been recently reported for Pt with 
conclusive results: The presence of H2 molecules alters 
the conductance histograms of clean nanocontacts^. It 
is thus important to determine the extent of this in- 
fluence when interpreting conductance histograms under 
low-vacumm conditions or in the presence of gases. We 
considered here two electrodes in the form of pyramids 
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FIG. 5: Transmission versus energy for V=0 for the AU-H2- 
Au bridge shown in the inset. 



in the (001) direction with a H2 molecule anchored be- 
tween them (see inset in Fig. |5J. The orientation of the 
molecule has been chosen on the basis of the report in 
Ref . 1291 for Pt, although other alternatives have been re- 
cently proposed^. We have checked that this orientation 
is stable for a range of distances d between outer planes 
which are fixed in the relaxation. Finally we have com- 
puted the conductance for d ~8 A(shown in Fig. [SJ. 
As one can see the conductance at the Fermi energy is 
greatly reduced from the quantum of conductance Go due 
to the presence of the H2 molecule. A detailed analysis 
of this problem will be published elsewhere^ 
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